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ABSTRACT 


The  trailing  vortices  generated  by  the  control  planes  of  sub- 
marines give  rise  to  surface  signatures  in  the  form  of  scars  and 
striations. 

Two  counter-rotating  vortices  were  generated  in  a  novel  experi- 
mental system  and  their  interaction  with  the  free  surface  was  investi- 
gated. In  addition,  the  governing  equations  have  been  solved  through 
the  use  of  the  boundary-element  method  for  a  representative  Froude 
number.  The  results  have  been  expressed  in  terms  of  the  depth  of 
submergence  of  the  vortices,  their  mutual  induction  velocity,  and  the 
initial  vortex  separation.  It  has  been  shown  that  the  free  surface 
begins  to  deform  when  the  vortices  are  at  a  distance  of  about  one  ini- 
tial vortex  separation  from  the  free  surface.  The  height  of  the  maxi- 
mum deformation  is  attained  at  a  normalized  time  of  about  0.1.  when 
the  vortices  are  at  a  distance  of  about  0.5  bo  from  the  free  surface. 
The  elevated  part  of  the  surface  is  bounded  by  two  scars,  whose 
moUon  is  slaved  to  that  of  the  vortices. 
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I.    INTRODUCTION 

Vortices  and  vortex  wakes  have  become  a  major  theme  of 
aerodynamics  research  since  the  advent  of  the  large  aircraft  and  the 
understanding  of  their  evolution  required  an  examination  of  many  of 
the  fundamental  problems  in  fluid  mechanics.  Much  of  the  progress 
made  during  the  past  two  decades  was  discussed  at  the  Symposium  on 
Aircraft  Wake  Turbulence  and  Its  Detection  [Ref.  1]  and  the  Aircraft 
Wake  Vortices  Conference  [Ref.  2).  Comprehensive  reviews  of  the 
entire  subject  have  been  given  by  Donaldson  and  Bilanin  [Ref.  3], 
Widnall  [Ref.  41,  and  Hallock  and  Eberle  [Ref.  5]. 

These  studies,  as  well  as  numerous  others  carried  out  since  1977, 
have  uncovered  a  number  of  complex  problems  which  must  be 
resolved  in  order  to  achieve  a  better  understanding  of  the  important 
features  of  trailing  vortices  in  homogeneous  and  stratified  media.  The 
principal  ones  are  as  follows  [Ref.  61: 

1.  Roll-up  process:  The  velocity  and  turbulence  distribution  at  any 
station  behind  the  wing  depend  on  the  wing  section,  wing-tip 
shape,  Reynolds  number,  wing  incidence,  and  the  distance  of 
the  station  from  the  wing  [Ref.  71.  The  distributions  of  the  initial 
velocity  and  turbulence,  which  influence  the  roll-up  and  the 
decay  process,  cannot  be  changed  independently.  For  example, 
a  change  in  the  tip  shape  changes  the  core  size,  as  well  as  the 
velocity  and  turbulence  distributions.  High  levels  of  turbulence 
result  in  an  increased  diffusion  of  vorticity,  which  in  turn 
increase  the  core  size. 

2.  Probe  sensitivity  of  the  vortices:  Flow  visualization  studies 
suggest  that  trailing  vortices  are  extremely  sensitive  to 
disturbances  created  by  even  very  small  probes  or  bubbles.   This 
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forces  one  to  use  non-intxusive  means  of  measurement  such  as  a 
Laser  Doppler  Velocimeter.  Even  then,  "vortex  wandering" 
[Ref.  8],  which  makes  the  vortices  appear  larger  than  normal  in 
time-averaged  velocity  measurements  (for  vortices  generated  by 
a  wing  in  a  wind  tunnel),  or  the  unsteady  nature  of  the  flow  (for 
vortices  generated  by  a  wing  in  a  tow  basin)  makes  the  mean 
velocity  profiles  in  the  vortices  difficult  to  determine. 

3.  Large-scale  instabilities:  The  vortices  are  seldom  observed  to 
decay  away  owing  to  viscous  and  turbulent  dissipation,  but  are 
almost  always  destroyed  by  either  mutual  induction  instability 
(Crow  instability  [Ref.  9])  and/or  vortex  breakdown.  The  Crow 
instability  grows  exponentially,  and  results  either  in  linking  of 
the  vortex  pair  into  a  series  of  crude  vortex  rings  or  in  a  highly 
disorganized  intermingling  of  the  vortices. 

Vortex  breakdown,  whose  mathematical  details  have  not  yet 
been  adequately  treated,  rearranges  the  vortex  structure  and 
increases  the  core  size,  turbulence,  and  energy  dissipation. 
Thus,  it  is  very  difficult  to  measure  accurately  the  trajectories  of 
the  three-dimensional  vortices  from  their  creation  to  their 
ultimate  demise. 

4.  Reynolds  number:  Even  the  highest  Reynolds  numbers,  based 
on  wing  chord,  reached  in  wind  tunnels  or  towing  basins,  are  an 
order  of  magnitude  lower  than  what  is  possible  for  an  aircraft. 
Thus,  the  scale  effects  are  not  easy  to  assess. 

5.  Ambient  conditions  such  as  turbulence  and  stratification  play 
major  roles  in  the  evolution  of  vortices.  The  quantification  of 
these  effects  requires  numerical  analysis  and  extremely  careful 
experiments  [Ref.  lOJ. 

6.  Ground  or  free  surface  effects:  The  vortex  pair  may  move  toward 
a  rigid  boundary  at  which  the  no-slip  condition  must  be  satisfied 
or  toward  a  free  surface  at  which  the  zero-shear  condition  must 
be  satisfied.  In  either  case,  the  vortices  come  under  the 
influence  of  their  images  and  move  accordingly. 

The  phenomenon  is  further  complicated  by  several  additional 

facts.    When  the  vortices  are  propelled  toward  a  rigid  surface,  vorticity 

of  opposite  sign  is  generated  on  the  no-slip  boundary  and  swept 

toward  the  vortex  pair.    The  total  vorticity  diminishes  very  quickly  as 

vorticity  from  the  two  regions  diffuses,  the  wall  region  serving  as  a 
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strong  sink  for  the  vorticity  associated  with  the  original  vortex 
[Ref.  11].  The  development  of  a  boundary  layer  along  the  rigid  wall 
may  give  rise  to  flow  separation  for  sufficiently  high  Reynolds  numbers. 
With  or  without  such  a  separation,  howfever,  the  center  of  the  vortex 
pair  eventually  moves  away,  or  "rebounds,"  firom  the  wall  IRefs.  11- 
13].        .  - 

For  the  case  of  a  zero-stress  boundary,  the  free  surface  still  acts  as 
a  vorticity  sink,  but  this  is  relatively  weak  due  to  the  absence  of 
intense  oppositely  signed  vorticity.  Thus,  in  the  absence  of  other 
impending  phenomena,  one  expects  a  mild  interaction  between  the 
vortices  and  the  free  surface  and  a  small  rebound  of  the  vortex  pair 
from  the  free  surface.  However,  the  ability  of  the  free  surface  to 
deform  under  the  influence  of  strain  fields  leads  to  a  strong  inter- 
action between  the  vortices  and  the  free  surface. 

It  is  evident  from  the  foregoing  that  the  motion  and  the  life-span 
of  trailing  vortices  are  governed  by  a  number  of  nonlinearly  dependent 
complex  phenomena.  A  number  a  experimental  and  analytical  studies 
have  been  carried  out  at  the  Naval  Postgraduate  School  by  Sarpkaya 
and  his  students  [Refs.  6,  14-211  in  order  to  investigate  the  effects  of 
these  parameters  on  the  rise  and  demise  of  the  trailing  vortices  in 
homogeneous  and  density  stratified  media.  These  studies  have  clearly 
identified  the  various  demise  mechanisms  in  both  media  and  estab- 
lished basic  relationships  between  the  rise  of  vortices  and  the 
governing  parameters  in  a  finite  as  well  as  effectively  infinite  medium 
[Ref.  201,  free  from  ambient  turbulence. 
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The  present  investigation  is  a  continuation  and  refinement  of  the 
previous  studies.  The  intent  is  to  analyze  the  rise  and  demise  of  a 
voriiex  pair  in  a  medium  with  a  deformable  free  surface.  This  problem 
is  of  interest  both  from  the  standpoint  of  determining  the  interaction 
effect  of  the  free  surface  on  the  rise  and  demise  of  the  vortex  pair  and 
from  the  standpoint  of  predicting  the  resulting  free  surface  shape.- 
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II.    PHYSICAL  ANALYSIS 

A.      DIMENSIONAL  ANALYSIS 

The  dependent  parameter  of  major  importance  to  the  problem 
solution  is  the  instantaneous  position  of  the  vortex  pair  (x,  y).  It  "may 
be  expressed  as  a  function  of  the  following  parameters  [Ref.  22]: 

X  =  f(t,  Vo,  do,  po.  dp/dy,  v,  bo,  g.  rg,  e,  Ln)  (1) 

and 

y  =  f(t.  Vo.  do.  po.  dp/dy,  V,  bo.  g,  re,  e,  Lii)  (2) 

in  which  the  variable  definitions  are  as  follows: 

t  time 

Vq        initial  mutual  induction  velocity  of  the  vortices 

do         initial  depth  of  the  vortex  pair 

po         reference  density  of  the  medium 

dp/dy  linear  density  gradient 

V  kinematic  viscosity  of  the  medium 

bo         initial  separation  of  the  vortex  pair 

g  gravitational  acceleration 

re  effective  core  radius  of  the  vortex 

e  rate  of  decay  of  the  turbulent  energy  per  unit  mass 

Li  1  '  integral  scale  of  the  turbulent  field 

The  height  and  width  of  the  test  section  were  not  included  in  the 

foregoing  because  a  detailed  analysis,  based  on  ideal  vortices,   has 

shown   that   the   velocities   induced   by  the   bottom    or   sides   were 
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negligible.  Effects  of  surface  tension  on  the  instantaneous  position  of 
the  vortex  pair  are  deemed  negligible  and  thus  are  not  included  as  a 
parameter  in  equations  (1)  and  (2). 

A  dimensional  analysis  of  equations  (1)  and  (2)  yields: 

x/bo  =  f(Vot/bo.  do/bo.  Nobo/Vo.  Vo2/gbo,  Vobo/v.  re/bo.  8*.  Ln/bo)^  (3) 
and 

y/bo  =  f(Vot/bo,  do/bo.  Nobo/Vo,  Vo2/gbo.  Vobo/v,  re/bo,  £*,  Li i/bo)     (4) 
in  which 

No  =  f^T^I  (5) 


'°-',Po     dyj 


is  known  as  the  Brunt-Vaisala  Frequency. 

All  parameters  in  equations  (3)  and  (4)  may  be  changed  indepen- 
dently except  re/bo,  which  is  taken  as  nature  provides  it.  The  primary 
reason  for  this  is  that  a  century  of  theoretical  and  experimental  aero- 
dynamics research  has  been  incapable  of  describing  the  details  of  the 
structure  of  the  tip  vortex  to  be  used  as  the  initial  conditions  in  the 
viscous  solution.  It  is  surprising,  but  true,  that  until  recently  the 
importance  of  the  generating  surface  shape  and  its  influence  upon 
both  the  initial  tangential  velocity  profile  and  the  initial  turbulence  in 
the  vortex  core  had  not  been  fully  appreciated.  Here  the  said  influ- 
ence has  been  characterized  in  terms  of  an  effective  core  radius  with 
full  awareness  of  its  shortcomings. 
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B.     GOVERNING    DIFFERENTIAL    EQUATIONS    AND    BOUNDARY 
CONDITIONS 

The  generation  of  Internal  waves  and  the  rise  and  demise  of  a 
vortex  pair  in  a  stratified  medium  may  be  analyzed  through  the  use  of 
the  equations  of  motion  for  an  incompressible  fluid.  These  equations 
may  be  applied  for  both  laminar  and  turbulent  motions  provided  that  a 
suitable  turbulence  closure  model  is  used  and  the  usual  Boussinesq 
approximation  (gravitational  acceleration  is  much  larger  than  fluid 
acceleration)  is  adopted.  For  the  type  of  motion  considered  herein, 
the  Boussinesq  approximation  is  quite  valid  and  has  been  used  in  the 
investigation  of  all  types  of  internal  waves  in  stratified  fluids. 

For  a  two-dimensional  flow,  with  y  vertical  and  x  horizontal,  the 
Navier-Stokes  equations  of  motion  are: 

dll  dU  dU  I     dP         'o  ,ry^ 

at  +  ^a^-^^5^=-Fa^^^^^  ^^^ 

av     av     av        1  ap     _« 

at  +  ^a^^^  =  g-F¥^''^^  ^^^ 

Differentiating  equations  (6)  and  (7)  with  respect  to  y  and  x  respec- 
tively yields 

a^u  a^u   a^u  avau  a2u_  i  a^ap  i_a2p   ^2  ^  to^ 
ayat  "^  axay  "^  ^axay  "^  ay  ay  "^  ^^  -  p^" ay  ax  p  axay  "*"  ^ay  ^^  ^^  ^^^ 

a^v  auav  a2v  _a^  a2v  _  1  a£_ap_j_a2p_  _a_,^2  ^  (a\ 
Bxdt  ■*"  ax  ax  "^  ^a^"^  axay  "*"  ^^^a^  -  p^'ax  ay  p  axay  "*■  ^ax  ^^  ^^  ^^^ 

Subtracting  equation  (9)  from  equation  (8)  yields 
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Subtracting  equation  (9)  from  equation  (8)  yields 


fdu     dv\ 


at 


dy    dx) 


dx 


y[dy   axJJ  ■*•  ay  L  \ay   ax] 


^  1  5£.5P 

""  p^ay  ax 


1  apap     ^2 


/^au_av^ 
ay   ax 


(10) 


1  ap 

Solving  equation  (6)  for  —  ^  yields 


1  aP       „o       Du 


(11) 


where 


D( )    in.  in.  Kl 

Dt 


at  ■•"^^IF"*""  ay 


1  ap 

Solving  equation  (7)  for  —  -^  yields 


1  ap  „o      Dv 


P  ^ 


Dt 


(12) 


Substituting  the  results  of  equations  (11)  and  (12)  into  equation  (10) 
and  defining  vorticity  as 


fau_av> 
^  =  lay   ax, 


(13) 
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results  in  the  following  expression: 


3^   acu^    a(vQ       2     £3? 


at""    dx 


By 


p  dx      LP  ^ 


l^M^^.m 


Dt 


pdyl  Dt 


P 


(14) 


When  the  gravitational  acceleration  is  several  orders  of  magnitude 
larger  than  the  fluid  accelerations,  the  Boussinesq  approximation  is 
appropriately  invoked  and  the  terms  in  brackets  in  equation  (14)  may- 
be neglected. 

In  order  to  deal  with  the  case  of  a  density  stratified  medium,  the 
density  is  defined  as 


P  =  Po+  p(y)  +  p'(x,y.t) 


(15) 


where  po  is  the  reference  density,  p  (y)  is  the  initial  variation  of  density 
with  y,  and  p'(x,y,t)  is  the  fluctuating  part  of  the  density  with  time. 
Then 

dx      dx 


and  equation  (14),  neglecting  the  bracketed  terms,  becomes 


3^   a(ua   a(vQ  x^e: 

At  +      ;^Y      +      Atr      -  vv    (,  -. 


dt^    dx 


dy 


Po  ax 


(16) 
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The  last  term  above  represents  the  effect  of  the  density  gradient  and 
gives  rise  to  oppositely  signed  vorticity  in  a  nonhomogeneous  fluid. 
The  diffusion  of  density  is  given  by 

Substituting  equation  (15)  for  density  in  the  above  equation  jrields 

The  equation  of  continuity  is 

Adding  equation  (19)  to  the  left  side  of  equation  (18)  and  simplifying 
results  in 

^*¥-^=-vt-VV.v0  (20, 

Equations  (16)  and  (20)  are  thus  the  governing  differential  equa- 
tions for  the  motion  of  a  vortex  pair  in  a  density  stratified  fluid. 

For  the  development  of  boundary  conditions,  it  is  assumed  at  the 
outset  that  the  influence  of  the  sidewalls  and  bottom  of  the  test 
section  on  the  rise  and  demise  of  the  vortex  pair  is  negligible.  The 
validity  of  this  assumption,  when  used  in  conjunction  with  numerical 
computations,  has  been  demonstrated  by  previous  investigators  [Refs. 
18-201.    The  fluid  domain  can  thus  be  considered  to  be  bounded  only 
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by  a  free  surface  which  can  be  described  as  a  function  of  x  and  time  as 
follows: 

y  =  ii(x.t) 
where 

Ti(x,o)  =  0 

describes  the  initial  undisturbed  location  of  the  free  surface. 

The  free  surface  requires  both  a  kinematic  and  a  dynamic  bound- 
ary condition  as  described  by  Sarpkaya  and  Issacson  IRef.  221.  The 
kinematic  condition  states  that  any  particle  which  lies  on  the  free 
surface  at  any  instant  will  never  leave  it.  This  leads  to 

Dt  =  at-*-^a]F  =  ^  ^^  y  =  ^  ^21) 

The^  dynamic  free  surface  condition  requires  that  the  pressure 
difference  across  the  interface  results  in  a  force  normal  to  the  bound- 
ary which  is  due  wholly  to  surface  tension.  This  condition  takes  the 
form 


P  =  Pa  +  a 


(—    —] 


(22) 


where  a  is  the  surface  tension,  b~  and  b"  are  the  radii  of  curvature  of 

Ri  R2 

the  free  surface  in  any  two  orthogonal  directions,  and  P  -  Pa  is  the 
pressure  difference  across  the  interface. 
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As  noted  previously,  the  effects  of  surface  tension  can  be  assumed 
to  be  negligible.  This  observation  has  been  verified  experimentally  by 
Gray  IRef.  17).  Thus,  if  the  flow  is  considered  inviscid,  the  pressure  P 
within  the  fluid  can  be  described  by  the  unsteady  Bernoulli  equation  as 
follows: 

where  0  is  the  velocity  potential  such  that 

"  =  ar  ^  =  -5r  (24) 

and  q2  =  u2  +  v2.   F(t)  is  an  arbitrary  function  of  time  only. 

When  the  pressure  just  outside  the  liquid  is  constant  (i.e.,  atmo- 
spheric), the  free  surface  condition  reduces  to 

90     q2 

30 
where  F(t)  has  been  included  as  part  of  -rr-. 

C.     DIMENSIONLESS   PARAMETERS 

It  is  convenient  at  this  point  to  cast  the  governing  equations,  (15) 
and  (20),  in  nondimensional  forms,  scaling  each  variable  by  a  quantity 
characteristic  of  its  expected  magnitude.  Two  possible  time  scales 
exist.  The  dynamic  time  scale  utilizes  the  time  a  characteristic  length 
would  be  traversed  by  a  fluid  particle  traveling  at  the  characteristic 
velocity.    The  buoyant  time  scale  is  based  on  the  natural  buoyancy 
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frequency  of  the  stratified  flow,  i.e.,  the  Brunt-Vaisala  frequency  Nq 
defined  by  equation  (5).    Each  of  these  scales  gives  a  slightly  different 
form  of  the  normalized  governing  equations. 
1.      Dynamic  Scale 

Introducing    Uc  and  Lc  as  the  characteristic  velocity  and 
length,  one  has  the  following  nondimensionalized  quantities: 


Cm  =  CLc/Uc  tm  =  Uct/Lc 


Um  =  u/Uc 


Vm  =  v/Uc 


Xm  —  x/Lc 


ym  =  y/Lc 


Pm  =  p/po 


(26) 


Substituting  the  above  into  equation  (16)  yields 

Uc    d     rCmUc^       1      d     (      ,,    CmUc^      L.J^(      tt 
Lcatml,   Lc    J  "■  Lc  axm  1,"°^^^     Lc    J  ""  Lc  aym  t"""^^' 

1    „2  /"CmUc^      g     a    ,     .  , 


Lc 


(27) 


Simplifying 


Uc^aCm      Uc^   3(UmCm)      Uc^   a(VmCm)      Uc     „2  r         £-^ 
L?'atm''L?'      axm      "^  L?      aym      "  L? '' ^m  ^m  +  j^  9x^ 


which  becomes 


aCm      a(UmCm)       a(VniCm) 


at 


m 


ax 


m 


-iiiZiiiL__JL_\7  2  r      ,  g^c  ap'm 
aym      "  UcLc  ^m  «>m  +  ujl  3xm 


(28) 


or 


aCm      a(UniCm)      a(vniCm) 


at 


m 


ax 


m 


ay: 


m 


Jl.t72  r      ^    1    ap'm 
R^V^Qm  +  F^axn, 


(29) 


25 


where 


Re  =  UcLc/v        (characteristic  Reynolds  number) 


Fv  =  Us/VgLc      (characteristic  Froude  number) 


Sim^ilarly,  the  density  diffusion  equation  may  be  expressed  in 
nondimensional  terms  by  making  equivalent  substitutions  into  equa- 
tion (20)  as  follows: 


do'    3(up')    a(vp')        dQ_     „2  .     ^_!5. 


at  ^   dx 


dy 


This  becomes 

L^at^  ^P°P'^  "^  Ca^  (UcUmPop'm)  +  E^^  (vmUcpoP'm) 


UcVm     a 


V  a2  (pop'm) 


LT^  (PoP'm)  +  iT  V^  (poP'nJ  +  T?^—^^ 


Which  upon  simplification  reduces  to 


UcPo  ^P'm      UcPo  a(UmPm)       UcPo  a(VmP'm) 
Lc     atm         Lc         axm  Lc        ayi 


^m 


vpoa^pm 


UcPo  apm  VPO  ^o        ,  -rw-r- 

Lc     ""^aymlT^niPm  +  l^-^ 


Multiplying  both  sides  by  Lc/poUc  yields 


apm      a(umP'm)        a(vmP'm) 
+  — ^n + 


at 


m 


ax 


m 


ay] 


m 


apm 


aym      UcLc 


V2pm  + 


a2pm^ 
"5^ 


(30) 
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Incorporating   t±ie   characteristic   Reynolds   number,    this   equation 
becomes 


3p'm      3(UmP'm)        3(vmp'm) 


ap 


m 


at 


m 


ax 


m 


ayi 


m 


-     """^  aym  ""  Re 


V  2  p  m  + 


a^p   ^ 


'm 


^, 


(31) 


The  Brunt-Vaisala  frequency  given  by  equation  (5)  as 


(    ^     ->-->l/2 


No  = 


-g  ae_ 

Po  ^y 


may  be  written  as 


o  "  poLc       aym 


which  becomes 


N^Lc       ap 


m 


g 


3y] 


m 


To  further  simplify  equation  (31),  the  following  definitions  are  made: 


N?  =  g/Lc 


(characteristic  Brunt-Vaisala  frequency 
squared) 


(33) 


and 


ap 


m 


5y] 


m 


(34) 


Substituting  equation  (34)  into  equation  (31),   the  nondimensional 
form  of  the  density  diffusion  equation  in  final  form  becomes 
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at 


m 


dx 


m 


ay. 


m 


^    \ 


•'^T —  +  — =v:: +  — ^r: =  ^'^m  +  s~  '  — 


'm-r  j^|V^p'in+  3y 


^ 


2 
my 


(35) 


2.     Buoyant  Scale 

Again  introducing  Uc  and  Lc  as  the  characteristic  velocity  and 
length,  the  following  nondimensional  quantities  are  defined: 


Vxn  =  u- 


_  x_ 


p-=^ 


In  this  case,  however,  time  is  nondimensionalized  using  the  square  of 
the  characteristic  Brunt-Vaisala  frequency.    Namely 


N^c 

=  g/Lc 

and  then 

tm 

=  tNc 

also 

._ 

Cm  =  C/Nc 

Substituting  the  above  into  the  governing  equation  given  by  equation 
(16)  yields 

a  1    a  1    a 

Nc  ^  (UNc)  +  L^a^  (umUcCmNc)  +1^^  (vmUcUNJ  = 


4vV2(UNe)+|:^3|^(poP'm) 
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simplifying 


^,2^U    UcNc  a    ,    .  ,    UcNc  a   ,    .  ; 


Nc     ^2r  g  apm 


which  becomes 


aCm        Uc    a(UmCm)        Uc    a(VmCm) 
atm  "*"  LcNc      axm      "^  LcNc      aym 


__v         2r  g     ap'm 


136) 


but 


Thus  equation  (36)  becomes 


aCm        Uc    ra(UmCm)      a(VmCm)^  V         ^c   „oy,         ^9m 


m 


ViL^l     3Xm      "      aym     j-Vii7l*"=    -msm-ax^ 


or 


at,   +^^ 


Lm 


aCUmCm).     a(Vm^m) 


ax 


m 


ay. 


m     J 


V|,Cm  + 


ap'i 


m 


Re     m'*m  -r  3x, 


m 


(37) 


where  Fy  and  Re  are  as  defined  above. 

The  density  diffusion  equation  is  similarly  nondimensional- 
ized  as 


321   a(u£i   a(ve2_     aj.      2       a£5 
at*    ax   *   av   -^  av  +  ^^  p +^  a7^ 


ay- 
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which  becomes 


The  above  equation  can  be  simplified  to  yield 

-,    9p'm      UcPo     d     ,         ,   .      UcPo     a     .        ,   , 

P°^^  at^'-TTa^^^^P^ ''  Tr^^'^^^P"^  = 

-UcPo        apm      PoV_2    ,        vpo  a^Pm 

"T^^-a^^T^^mPm  +  i^-^^ 

Introducing  Fv  and  Re  as  before  and  recalling  that  Nj  =  g-,  one  has 


ap' 


m 


at 


m 


+  Fv 


'a(Ump'm)      a(vmp'm) 
.     axm  aym 


=  FvVin^c^  + 


ay. 


m 


Re 


/ 


v^pv 


a2p  ^ 


'm 


(38) 


Now  as  developed  above  in  equation  (34) 


n2  = 


ap 


m 


ay: 


where  n^  = 


N2      N^Lc 


m 


nT  = 


g 


a  Pm 
Substituting  for  ^- —  in  equation  (38)  yields  the  nondimensionalized 

density  diffusion  equation 


ap', 


m 


at 


+  F^ 


m 


'a(Ump'm)       a(Vmp'm) 


ax 


m 


ay. 


m     J 


=   Fvn2vm  + 


Re 


T72    .        ^~Pm 


5y. 


m 


(39) 
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Equations  (37)  and  (39)  are  valid  when  Fv«l  and  buoyancy 
dominates  the  flow.  When  Fy  approaches  zero,  the  equations  that 
result  from  equations  (37)  and  (39)  describe  the  propagation  of  linear 
internal  waves.  The  buoyant  scaled  forms  of  the  governing  equations 
are  of  interest  to  the  present  investigation  because,  for  submerged 
bodies  of  naval  interest,  Fy  is  typically  about  0.001  and  thus  signifi- 
cantly less  than  one. 

Having  established  the  buoyant  scale  as  the  form  of  interest, 
the  boundary  conditions  can  also  be  nondimensionalized  using  the 
same  technique.    From  equation  (21), 

ari        art 

3^+U33^  =  v  aty  =  Ti 

which  becomes 

where 

Tim  =  Tl/Lc 

Simplifying  the  above  jrields 

^Jlm        Uc  ^Tlm        Uc 

atm  "^  NcLc  ^^axm  "  NcLc  ^"^ 

Introducing  Fy  as  before 

3^  +  Fy  Um  ^^  =  FyVm  at  ym  =  Tim  (40) 
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and 


0 
0m  =  AT  T  2        (nondimensional  velocity  potential) 

C       Q 


Equation  (25)  can  be  nondimensionalized  as  follows: 
30     q2 


^  +  :^+gn  =  o 


becomes 


Nc atm         +  "T^  •<•  gLcnm  =  0  (41) 


where 


qm  =  ^  +  v^ 


Further  simplifying  equation  (41)  yields 

a0m  -  U?   q2^     gn 
atxn  "•  N^  2    +  N^ 

but 


""  =  0  (42) 


N-^Lc=l  and  ^^  =  n 


Thus  equation  (42)  becomes 


a0m      ^2    ^m  ^        , 

•^:^+F5  — +  llm  =  0     at    ym=T|m 
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For  the  purposes  of  tJie  present  investigation,  the  flow  will  be 
assumed  to  be  inviscid,  i.e.,  the  viscous  difi'usion  will  be  ignored  to  a 
first  order  approximation.  The  nondimensionalized  forms  of  the 
governing  equations  and  boundary  conditions  are  summarized  in 
Figure  1. 
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m.    NUMERICAL  SOLUTION  TECHNIQUES 

A.     INTRODUCTION 

The  significant  difference  between  this  study  and  the  work  of 
previous^^  investigators  is  the  introduction  of  a  deformable  free  surface 
and  the  resulting  complex  interplay  between  kinetic  and  potential 
energies.  The  computational  domain  thus  has  an  unknown  boundary 
on  which  a  double  condition  must  be  imposed  as  previously  outlined. 
This  complexity,  as  well  as  several  other  specific  features  of  this 
problem,  directly  influences  the  ability  of  a  numerical  method  to  con- 
verge to  an  adequate  solution. 

The  governing  differential  equations  and  boundary  conditions  are 
both  nonlinear.  Sarpkaya  and  Issacson  [Ref.  23]  note  that  for  small 
amplitude  waves  (amplitude  «  wavelength),  the  boundary  conditions 
at  the  free  surface  may  be  linearized.  This  approach,  although  possibly 
valid  in  the  vicinity  of  surface  striations,  would  not  be  valid  for  the  scar 
front  where  observed  light  diffraction  patterns  [Ref.  171  indicate  very 
small  radii  of  curvature.  Additionally,  Haussling  and  Coleman  [Ref.  24] 
have  demonstrated  numerically  the  importance  of  nonlinear  terms  in 
solving  the  free  surface  problem  in  the  vicinity  of  the  generation 
source.   This  is  the  case  encountered  in  the  region  of  the  scar  front. 

The  rate  of  change  of  the  free  surface  deformations  also  tends  to 
be  slow.  The  scar  pattern  develops  as  the  vortex  pair  rises  to  its 
maximum  height  and  then  the  scar  is  trapped  by  and  slaved  to  the 
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vortex  beneath  the  surface.  The  scars  and  striations  are  thus  not  small 
amplitude  surface  waves  propagating  independently  from  their 
generator,  but  are  in  fact  local  disturbances  whose  generation,  growth, 
and  propagation  are  controlled  by  the  vortices.  This  aspect,  observed 
experimentally  by  Gray  [Ref.  17],  separates  this  problem  from  other 
research  in  the  field  of  deformable  free  surfaces.  Previous  investiga- 
tors, applying  numerical  methods  to  such  surfaces,  have  dealt  almost 
exclusively  with  surface  waves  generated  by  a  moving  body  (surfaced  or 
submerged)  and  propagating  away  independently. 

The  numerical  methods  reviewed  for  solution  of  this  problem  can 
be  broadly  categorized  into  Finite  Differences,  Finite  Elements, 
Boundary  Integral  Equations,  and  Hybrid  Methods  which  employ  com- 
binations of  the  others.  For  the  purposes  of  this  study,  such  catego- 
rization is  based  on  the  manner  in  which  the  governing  equations  are 
tackled  in  the  computational  domain.  It  is  noted  that  all  of  the 
methods'  reviewed  utilize  some  form  of  finite  difference  scheme  with 
respect  to  time.  The  applicability  of  each  method  and  the  associated 
advantages  and  disadvantages  are  discussed  and  summarized  at  the 
end  of  this  chapter. 

B.     FINITE  DIFFERENCE  TECHNIQUES 

The  use  of  Finite  Difference  techniques  in  solving  partial  differen- 
tial equations  has  been  well  established  and  successfully  implemented 
by  several  investigators.  For  the  case  of  a  nondeforming  free  surface, 
the  governing  equations  are  solved  in  the  computational  domain  using 
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at  ^  dx  ^  ay  -  ^^  ^    Po  ax 


and 


ap'    a(up)    a(vp)      vap     „2  .      ^^p 


and  the  Biot-Savart  integrals 


(30) 


(y'-y)C(x'.y')dx'dy' 


A 


,      ,        f  (X  -X')  C  (x\  y')  dx'dy'  .... 

v(x,y)  =    J  ^^^;^ ^  (32) 

A 

where 

r2  =  (x  -  x')  +  y  -  y') 

x',  y*:    position  of  the  vorticity 

X,  y:      point  where  u  and  v  are  to  be  calculated 

In  nondimensional  form,  the  buoyantly  scaled  forms  of  equations  (31) 
and  (32)  become 


36 


r  uj> 


Um  (Xm.ym)  = 


(y'm  -ym)  Cm  (x'm.  y'm)  dxmdy'm  LcNc 


J 
A 


2^ 


Ur 


Vm  (Xm.ym)  = 


'  (Xjn-x'm)  Cm  (x  m.  y'm)  dx'mdy'm  LcNc 


27cr 


m 


Uc 


Upon  simplification,  the  equations  become 


Fv  Um  (Xm.ym)  - 
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where  parameters  are  nondimensionalized  as  before. 

Thegoveming  equations  can  then  also  be  rewritten  as 
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where  FyUm  and  FyVm  are  calculated  directly  using  the  nondimension- 
alized forms  of  the  Biot-Savart  Integrals  equations  (33)  and  (34). 
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The  treatment  of  the  governing  differential  equations  in  a  stream 
function  formulation  is  thus  relatively  straightforward.  The  complica- 
tion is  associated  with  the  determination  of  the  free  surface  location  in 
conjunction  with  the  solution  of  the  governing  equations,  all  preferably 
simultaneously,  as  noted  by  Yeung  [Ref.  25]. 

Finite  Difference  methods  are  most  suitable,  or  at  least  simplest 
to  implement,  for  a  boundary  geometry  that  is  rectilinear.  The 
deformable  free  surface  requires  the  use  of  "irregular  stars"  for  differ- 
ence schemes  near  the  fluid  boundary.  These  "irregular  stars"  are  of 
an  inferior  accuracy  when  compared  to  the  remainder  of  the  grid. 
Thus,  either  the  mesh  must  be  refined  or  the  difference  formula 
changed  to  a  higher  order  if  accuracy  is  to  be  maintained  near  the  free 
surface  and  especially  in  the  vicinity  of  the  scars.  Yeung  [Ref.  25] 
notes  that,  with  the  location  of  the  free  surface  unknown  a  priori,  we 
have  the  unfortunate  situation  that  the  regions  that  demand  the  great- 
est accuracy  are  precisely  those  where  it  is  hardest  to  achieve.  Also, 
since  the  shape  of  the  free  surface  affects  the  migration  of  the  vortex 
pair,  a  loss  of  accuracy  in  determining  the  free  surface  will  be 
reflected  in  the  determination  of  the  location  of  the  vortex  pair.  This 
complication  is  one  not  included  in  the  investigations  reviewed  below 
where  the  generating  body  moves  independently  of  the  free  surface 
location. 

Conformal  transformations  can  simplify  one  aspect  of  the  problem 
by  transforming  the  real  geometry  with  a  deformed  free  surface  to  a 
rectilinear  mesh.    The  boundary  conditions  are  thus  simplified  but  the 
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resulting  field  equations  are  increased  in  complexity.  Haussling  and 
Coleman  have  successfully  utilized  transformation  functions  of  this 
form  to  generate  a  time  independent  computational  region  in  which 
nonlinear  free  surface  boundary  conditions  are  applied. 

The  Finite  Difference  method  primarily  benefits  from  the  direct 
solution,  of  the  governing  differential  equations  as  opposed  to-  the 
development  of  intermediate  integral  relations  required  in  the  inte- 
gral boundary  techniques.  The  cost  of  such  an  approach  is  in  the 
computational  intensity  required  to  iteratively  solve  the  finite  differ- 
ence forms  over  the  computational  domain.  Haussling  and  Coleman 
have  demonstrated  the  use  of  a  successive  over  relaxation  technique  to 
solve  steady  state  nonlinear  free  surface  boundary  condition  problems 
as  discussed  above. 

A  modification  of  successive  over  relaxation,  and  one  offering 
increased  rates  of  convergence,  has  been  demonstrated  by  Brandt. 
Dendy,  and  Ruppel  [Ref.  26].  Their  technique  utilizes  a  multigrid 
solution  which  solves  for  low-frequency  components  of  error  on  a 
course  grid  where  the  calculation  is  relatively  inexpensive,  and  high- 
frequency  components  of  error  on  a  fine  grid  where  successive  over- 
relaxation  is  efficient. 

Theodussiou  and  Sousa  [Ref.  27]  also  utilized  a  modified  grid  sys- 
tem to  speed  convergence  by  staggering  their  grid  structure  such  that 
pressure  was  defined  at  the  center  of  the  discretized  control  domain 
and  velocity  components  were  defined  at  the  center  of  the  control 
domain  faces.    This  arrangement  has  the  convenient  feature  that  the 
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velocity  components  are  stored  at  just  the  points  at  which  they  are 
required  for  the  calculation  of  their  advective  contributions.  The 
pressure  gradients  can  be  represented  by  central  differences  without 
inducing  non-physical  oscillations  in  the  pressure  distribution.  Note, 
however,  that  both  of  these  multigrid  techniques  add  to  the  complex- 
ity of  establishing  a  mesh  which  moves  with  the  deforming  free  sur- 
face. Both  sets  of  authors  suggest,  however,  that  this  added  complex- 
ity is  outweighed  by  the  savings  in  ease  of  convergence. 

Ohring  [Ref.  281  and  Ohring  and  Telste  [Ref.  291  have  also  partially 
circumvented  the  computational  intensity  of  an  iterative  technique  by 
directly  solving  the  finite  difference  equations  resulting  from  the 
Laplace  Equation.  The  technique  employed  utilizes  a  fourth-order 
solver  to  diagonally  decompose  the  resulting  coefficient  matrix.  Taylor 
series  approximations  are  also  suggested  for  application  to  the  nonlin- 
ear free  surface  boundary  conditions.  It  is  suspected,  however,  that 
the  computational  benefits  derived  from  such  a  technique  will  be  lost 
when  the  coefficient  matrix  becomes  nonlinear  as  would  result  from 
the  governing  equations  in  this  problem. 

For  completeness,  as  noted  previously,  finite  differencing  in  time 
is  required  for  all  the  numerical  methods  reviewed.  Sarpkaya  and  his 
students  [Refs.  17-201  have  successfully  employed  an  upwind-differ- 
encing scheme  in  time  and  verified  it  with  experimental  results. 
Yeung  notes  that  the  free  surface  conditions,  being  first  order  in  time, 
can  be  used  to  advance  the  solution  of  the  elevation  and  velocity 
potential  on  the  free  boundary.    However,  the  difference  form  utilized 
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can  dramatically  affect  the  stability  and  thus  a  modified  Euler  method 
is  suggested  since  it  is  unconditionally  stable.  Haussling  and  Coleman 
successfully  utilized  this  technique  for  advancing  their  solution  in 
time. 

C.  FINITE  ELEMENT  METHODS 

The  Finite  Element  and  Finite  Difference  methods  share  the 
common  feature  that  both  attempt  to  solve  the  governing  differential 
equations  directly,  differing  only  in  methodology.  The  Finite  Element 
Method  is  one  based  on  the  method  of  weighted  residuals.  The  usual 
procedure  consists  of  first  subdividing  the  domain  of  interest  into  a 
mesh  of  finite-sized  subregions,  within  each  of  which  the  solution  is 
represented  by  some  convenient  choice  of  trial  functions,  usually 
polynomials.  The  trial  functions  are  determined  by  substituting  them 
into  the  governing  equations  and  requiring  the  integrated  error  or 
residual  based  on  certain  weighing  functions  to  vanish.  An  integration 
by  parts  is  normally  preformed  to  reduce  the  interelement  continuity 
requirements  of  the  trial  functions  and  to  incorporate  the 
nonhomogeneous  boundary  conditions.  The  weighing  functions,  also 
known  as  test  functions,  can  be  chosen  in  a  variety  of  ways.  A  "weak 
formulation,"  such  as  the  Galerkin  Method,  makes  the  space  of  the 
test  functions  identical  to  that  of  the  trial  functions.  In  contrast,  a 
"strong  formulation"  is  one  based  on  the  existence  of  a  variational 
principle  where  a  functional  is  made  stationary.  Specifics  of  finite 
element  techniques  may  be  found  in  Zienkiewicz  [Ref.  30]  and  Dhalt 
andTouzot  [Ref.  311. 
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The  nonlinear  aspects  of  this  problem  affect  the  Finite  Element 
Method  in  much  the  same  way  as  the  Finite  Difference  Method.  The 
coefficient  matrix  of  a  representation  such  as 

[K]  [U]  =  IB] 
where 

[K]  =  coefficient  matrix  .  ~ 

[U]  =  matrix  of  unknowns 

[B]  =  forcing  function  matrix 
is  neither  symmetric  nor  independent  of  [U],  as  in  the  case  of  a  linear 
problem.  An  iterative  technique  is  thus  required  to  solve  this  problem 
at  a  given  instance  in  time.  As  a  result,  much  of  the  advantage  of 
reduced  storage  and  computation  time  inherent  to  finite  element  for- 
mulations is  lost.  The  choice  of  iterative  techniques  does  not  differ 
from  that  available  to  Finite  Difference  Methods  and,  thus,  once  the 
computational  domain  is  discretized,  no  significant  difference  in 
problem- solution  is  involved. 

The  discretization  of  the  computational  domain  is  an  area  where 
the  Finite  Element  Method  does  have  distinct  advantages.  The  intro- 
duction of  curvilinear  or  isoparametric  elements  of  higher  order 
allows  one  to  cope  with  any  arbitrary  boundary  geometry  with  little 
loss  of  accuracy.  Yeung  notes  that,  although  this  particular  advantage 
is  reduced  when  Finite  Difference  Methods  are  used  with  conformal 
transformations.  Finite  Element  Methods  still  retain  superiority  in 
flexibility,  particularly  in  the  case  where  varying  size  and  shape 
elements   are   introduced   to   overcome   local   irregularities.      Such 
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irregularities  are  exactly  the  problem  involved  in  the  vicinity  of  the 
deformable  free  surface.  Curvilinear  elements  could  be  used  to  fit  the 
free  surface  with  finer  elements  incorporated  in  the  vicinity  of  the 
scar  front.  This  could  be  accomplished  with  little  loss  of  accuracy  and 
less  complexity  than  that  which  results  from  the  use  of  "irregular 
stars"  in  the  Finite  Difference  Method.  Note,  however,  that,  since  the 
scar  front  moves  with  time,  an  adaptive  mesh  refinement  technique 
would  be  required  to  appropriately  place  the  finer  elements  in  the 
vicinity  of  the  scar.  Sarpkaya  and  Hiriart  IRef.  321  successfully  utilized 
varying  size  elements  in  the  vicinity  of  the  free  surface  in  conjunction 
with  their  moving  net  computation.  Although  the  free  surface  in  their 
case  differed  from  that  involved  in  this  problem,  the  basic  concept  of  a 
"flexible  mesh"  remains  unchanged. 

Admittedly,  the  basic  problem  of  determining  the  location  of  the 
free  surface  in  conjunction  with  solving  the  governing  equations 
remains  one  of  trial  and  error,  regardless  of  the  type  of  discretization 
employed.  Larock  and  Taylor  [Ref.  33]  adjusted  the  free  surface  loca- 
tion to  achieve  tangency  to  the  surface  velocities  calculated  based  on 
an  assumed  free  surface  position.  The  pressure  boundary  condition 
was  then  used  as  a  check  of  this  resulting  location.  Larock  and  Taylor 
note,  however,  that  such  a  technique  will  not  work  well  where  high 
curvature  or  Froude  numbers  (Fy)  less  than  one  are  involved,  as  is  the 
case  in  the  present  investigation.  As  an  alternative,  Sarpkaya  and 
Hiriart  adjusted  their  free  surface  location  to  satisfy  both  the  pressure 
and  velocity  boundary  conditions  simultaneously. 
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The  Finite  Element  Method,  when  formulated  using  the  varia- 
tional techniques,  also  benefits  from  the  fact  that  the  boundary  condi- 
tions may  be  included  as  an  integral  pari:  of  the  functional  rather  than 
dealt  with  separately.  Bai  and  Yeung  [Ref.  341  and  others  have  suc- 
cessfully employed  variational  techniques  to  directly  incorporate 
boundary  conditions  in  the  solution  of  linear  water  wave  problems,- 

D.     BOUNDARY  INTEGRAL  EQUATION  METHODS 

The  term  "Boundary  Integral  Equation  Method"  can  be  applied  to 
a  large  group  of  numerical  techniques  that  include  Green's  Function 
formulations.  Spectral  Methods,  and  a  boundary  element  application  of 
the  Finite  Element  Method.  In  all  cases,  the  solution  approach  differs 
dramatically  from  that  of  Finite  Differences  and  Finite  Elements  in 
that  the  governing  equations  are  not  attacked  directly.  Instead,  the 
problem  is  solved  by  satisfying  boundary  conditions  on  a  discretized 
boundary  and  thus  reducing  the  spatial  dimensions  of  the  problem  by 
one.  Physical  quantities  such  as  wave  height  and  fluid  pressures  are 
required  and  solved  for  only  on  the  boundaries.  Interior  data,  although 
available  based  on  the  boundary  solution,  is  not  specifically  required. 
Boundary  Integral  Equation  Methods  thus  have  the  distinct  advantage 
that  only  the  physical  quantities  specifically  desired  on  the  free  sur- 
face are  required  for  the  solution. 

Basically,  there  are  two  approaches  to  boundary  integral  formula- 
tions. Either  an  approximating  function  is  chosen  on  the  boundary 
that  satisfies  the  governing  equations  in  the  domain  and  approximates 
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the  boundary  conditions,  or  vice  versa.  In  both  cases,  the  solution 
technique  can  be  further  divided  into  indirect  and  direct  methods. 

In  indirect  approaches,  the  fundamental  solution  is  approximated 
on  the  boundary  by  a  function  with  unknown  coefficients.  These  coef- 
ficients are  found  by  satisfying  the  boundary  conditions.  Distributed 
singularity  methods,  such  as  simple-source  distributions,  are  an  exam- 
ple of  this  approach.  In  this  case,  the  fundamental  solution  is  repre- 
sented on  a  discretized  boundary  by  distributed  simple  sources  of 
unknown  strength  at  known  locations.  The  strength  of  each  source  is 
then  determined  based  on  the  solution  of  the  boundary  conditions. 

Direct  methods  find  the  fundamental  solution  through  the  use  of 
Green's  Function  formulations,  which  directly  incorporate  the 
governing  equations.  The  resulting  Green's  Function  integrals  are 
typically  solved  using  quadrature  or  point  kernel  techniques.  The 
main  disadvantage  in  this  approach  is  that,  for  the  governing  equations 
involved  in  this  problem,  there  is  no  straightforward  Green's  Function 
formulation  that  leads  to  fundamental  solution. 

The  main  advantage  of  the  Boundary  Integral  Equation  Method  is, 
then,  the  ability  to  directly  discretize  the  firee  surface  without  a  loss  of 
accuracy.  Information  is  thus  obtained  exactly  where  required.  How- 
ever, Brebbia  [Ref.  35]  notes  that  this  is  not  without  the  sacrifice  of  a 
symmetric  coefficient  matrix  such  as  that  which  is  common  to  the 
Finite  Element  Method.  Yeung  notes  that,  in  general,  the  great 
reduction  in  the  size  of  the  matrix  outweighs  the  added  complexity  of 
solving  a  nonsymmetric  system  of  equations. 


45 


The  disadvantages  of  Boundary  Integral  Equation  Methods  are 
directly  associated  with  errors  that  are  the  result  of  discretization. 
This  discretization  plays  such  a  large  role  in  the  formulation  of  an 
efficient  boundary  integral  equation  that  Brebbia  refers  to  these  tech- 
niques as  "Boundary  Element  Methods."  In  all  formulation  tech- 
niques, elements  of  one  type  or  another  are  formed  over  whi^h  a 
fundamental  solution  must  be  approximated  either  by  distributed  sin- 
gularities. Green's  Function  integrals,  or  trial  functions  for  finite 
elements. 

Discretization  errors  result  both  from  collocation  errors  and  geo- 
metric surface  errors.  Collocation  errors  are  primarily  the  result  of 
leakage,  as  noted  by  Hunt  [Ref.  36].  Since  the  boundary  conditions  are 
satisfied  exactly  only  at  discrete  points,  the  remainder  of  the  free  sur- 
face is  "porous"  by  comparison.  Leakage  errors  are  reduced  by 
increasing  the  number  of  free  surface  elements  and  thus  the  number 
of  discretization  points.  This,  of  course,  is  done  at  the  expense  of 
computational  intensity.  Geometric  surface  errors  result  from  the 
approximation  of  a  curved  free  surface  by  linear  elements.  This  prob- 
lem becomes  particularly  noticeable  in  areas  of  high  curvature,  such  as 
scar  fronts.  Higgins  and  Cokelet  [Ref.  37]  noted  that  this  problem  can 
be  partially  circumvented  by  using  a  Lagrangian  description  of 
"marked"  particles  on  the  free  surface.  These  particles  tended  to 
concentrate  in  the  regions  of  highest  curvature,  thus  giving  improved 
accuracy  exactly  where  needed. 
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Finite  Element  Method  formulations  using  boundary  elements 
only  can  also  circumvent  some  of  the  errors  associated  with  the  collo- 
cation and  geometric  surface  by  using  curvilinear  or  higher  order 
parametric  one-dimensional  elements.  The  problem  then  is  one  of 
finding  suitable  fundamental  solutions  that  satisfy  both  the  governing 
equations  and  continuity  requirements  between  elements. 

Finally,  Boundary  Integral  Equation  Methods  are  complicated  by 
the  necessity  to  include  nonlinear  terms  in  the  boundary  condition  for 
the  free  surface.  As  previously  noted,  this  is  necessary  in  order  to 
maintain  accuracy  in  the  vicinity  of  the  scar  front  where  very  small 
radii  of  curvature  occur.  Several  researchers  have  successfully  applied 
various  forms  of  the  Boundary  Integral  Equation  Method  to  Laplace 
equations  with  linearized  free  surface  conditions,  but  only  a  few  have 
incorporated  a  nonlinear  condition.  Faltinsen  [Ref.  381  incorporated 
nonlinear  conditions  by  using  the  properties  of  the  fluid  particles  on 
the  free  surface  at  one  instance  in  time  to  establish  a  new  free  surface 
location  and  step  the  solution  forward  in  time. 

E.     HYBRID  METHODS 

The  use  of  a  different  technique  in  different  portions  of  the  com- 
putational domain  results  in  a  hybrid  numerical  formulation.  This  type 
of  formulation  attempts  to  maximize  the  benefits  of  any  one  particular 
method  by  employing  it  only  in  regions  where  its  accuracy  remains 
high.  A  secondary  objective  is  to  reduce  the  computational  intensity  of 
the  overall  routine  by  using  a  coarser,  less  time-intensive  technique  in 
areas  where  accuracy  is  not  of  particular  concern.    This  is  exactly  the 
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case  in  the  regions  far  removed  from  the  free  surface.  The  usual 
approach  is  to  take  advantage  of  the  availability  of  analytical  solutions 
in  regions  where  the  flow  geometry  is  relatively  simple.  This 
approach  allows  the  number  of  mesh  points  to  be  reduced  with  a 
resulting  decrease  in  required  storage  and  computational  intensity.  A 
further  advantage  is  that  the  analytical  solutions  can  be  chosen  to 
permit  a  simple  solution  of  radiation  tjrpe  boundary  conditions,  as 
noted  by  Yeung  [Ref.  25]. 

The  disadvantages  of  Hybrid  Methods  are  associated  with  the 
necessity  to  properly  match  the  solutions  of  different  subdomains  in 
the  overall  computational  domain.  This  matching  may  result  in 
numerical  perturbations  to  the  solution  technique  if  not  properly 
employed.  Additionally,  the  advantage  of  reduced  total  grid  points  is 
somewhat  offset  by  the  added  computations  required  to  match 
solutions  at  the  common  boundaries. 

Yeung  notes  that  "the  more  successful  Hybrid  methods  have  so  far 
been  restricted  to  linearized  problems  where  analytical  solutions  in 
the  exterior  regions  could  be  obtained  without  too  much  difficulty.  In 
particular,  treatment  of  steady  flows  in  a  uniform  stream  or  time-har- 
monic flows  with  linearized  boundary  conditions  have  been  quite  well 
established."  [Ref.  25] 

Bai  and  Yeung  [Ref.  34]  utilized  a  finite  element  grid  in  the  vicinity 
of  the  free  surface  and  an  eigenfunction  solution  in  the  exterior  region 
to  reduce  the  computational  intensity  of  their  routine.  Chang  and  Pien 
[Ref.  39]  used  a  source  distribution  to  formulate  a  Boundary  Integral 


48 


Equation  MetJiod  near  the  free  surface  and  a  finite  difference  routine 
in  the  exterior  regions  to  calculate  the  hydrodynamic  forces  on  a  body 
moving  beneath  a  free  surface.  It  should  be  noted,  however,  that  both 
sets  of  investigators  used  a  Laplace  formulation  with  linearized  bound- 
ary conditions  to  obtain  their  results. 

For, completeness,  it  is  noted  that  Hybrid  methods  could  be  con- 
sidered to  include  various  numerical  techniques  that  employ  the  same 
method  throughout  the  computational  domain  but  in  different  man- 
ners in  each  subdomain.  Such  is  the  case,  for  example,  when  a  finite 
element  method  is  employed  with  a  varying  grid  and /or  element  type 
in  various  regions  of  the  computational  domain. 

F.      NUMERICAL  METHODS.   CONCLUSIONS.  AND  SELECTION 

The  sjniopsis  of  advantages  and  disadvantages  for  each  numerical 
method  listed  in  Table  1  provides  a  good  source  of  information  for 
making  a  wise  choice  of  formulation  to  be  used  in  this  problem.  The 
"best"  method  is  one  which  will  meet  the  goals  of  the  investigation 
while  maximizing  the  advantages  in  areas  of  particular  concern. 

In  this  case,  accuracy  in  the  vicinity  of  the  free  surface  is  of  par- 
ticular importance  since  ultimately  it  is  the  free  surface  shape  which 
is  desired.  This  accuracy  must  be  obtained  with  full  consideration  of 
the  necessity  to  include  nonlinear  boundary  conditions  while  mini- 
mizing the  computational  intensity  of  an  iterative  process.  Tuck,  in  a 
paper  by  Bai  and  Yeung,  observed  that  "to  a  certain  extent  the  *best' 
method  will  always  be  that  which  appeals  most  to  the  person  pro- 
gramming it,  and  hence  that  which  gives  him  the  greatest  chance  of 
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writing  a  successful  program  irrespective  of  efficiency.  There  is  no 
less  efficient  program  than  one  which  does  not  work  at  all."  [Ref.  34] 
Given  that  a  successful  program  can  be  developed  (a  pretty  big  given), 
it  is  then  necessary  to  maximize  both  efficiency  and  accuracy. 

With  full  consideration  of  the  concerns  given  above,  the  Boundary 
Integral' Equation  Method  surfaces  as  the  method  of  choice  for  the 
following  reasons: 

1.  Computational  intensity  is  minimized  by  reducing  the  spatial 
dimensions  by  one.  This  aspect  is  particularly  noticeable  when 
consideration  is  given  to  the  iterative  requirements  of  this 
problem.  This  substantial  improvement  directly  incorporates 
any  similar  advantage  that  can  be  gained  by  using  a  varying  Finite 
Element  Method  or  Hybrid  Method. 

2.  Accuracy  is  maintained  at  the  free  surface  by  incorporating  the 
best  aspects  of  Finite  Element  Methods  and  a  Lagrangian 
description  of  marked  particles.  This  direct  discretization  of 
the  free  surface  provides  both  accuracy  where  needed  most  and 
an  exact  solution  of  boundary  conditions  at  the  marked  points. 

3.  Storage  requirements  are  drastically  reduced  since  only  physical 
quantities  at  the  free  surface  are  required  or  needed. 

The  implementation  of  the  Boundary  Integral  Equation  Method 

through  the  use  of  distributed  vortices  will  be  outlined  in  detail  in  the 

following  section. 
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IV.    DISTRIBUTED  VORTEX  MODEL 

A.     INTRODUCTION 

Discrete  vortices,  with  or  without  a  core,  or  vortex  sheets  have 
been  used  as  boundary  elements  to  simulate  separated  and  unsepa- 
rated  flows.  The  method  consists  of  the  determination  of  the  appro- 
priate strength  of  the  vortices  at  each  time  interval  through  the  use  of 
the  governing  equations.  The  vortices  are  then  convected  to  their 
new  positions  and  the  process  is  repeated.  In  spite  of  its  simplicity, 
the  distributed  vortex  model  presents  several  difficulties,  all  of  which 
are  related  to  discretization  and  the  use  of  vortices.  The  evaluation  of 
the  governing  integral  equation  cannot  accurately  be  accomplished 
merely  by  applying  a  standard  integration  formula,  A  vortex  sheet  or  a 
string  of  point  vortices  is  unstable  to  small  sinusoidal  disturbances  of 
any  wavelength.  This  phenomenon  of  Helmholtz  instability  persists  in 
curved  nonuniform  vortex  sheets,  at  least  for  short  waves,  unless  the 
sheet  is  rapidly  stretching.  In  other  words,  the  round-off  and  trunca- 
tion errors  are  rapidly  amplified  to  cause  the  chaotic  motion  which 
often  ruins  practical  calculation.  If  it  is  granted  that  it  is  the  growth  of 
short  waves  which  can  ruin  calculations  with  vortex  sheets,  it  is  sensi- 
ble to  consider  ways  of  removing  the  instability.  This  is  because  the 
instability  is  introduced  by  the  step  of  replacing  a  shear  layer  of  small, 
but  finite,  thickness  by  a  vortex  sheet.  One  could  give  up  the  vortex 
sheet  approximation  and  return  to  the  computation  of  the  evolution  of 
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a  thin  layer.   This  is  wit±iout  doubt  the  most  satisfactory  procedure,  but 
it  involves  much  more  computation. 

An  alternative  approach  is  to  modify  the  governing  differential 
equation  to  allow  for  finite  thickness  but  the  resulting  equation,  while 
only  a  little  more  complicated  than  the  original  equation,  has  not 
proved  amenable  to  computation.  — 

Another  possibility  is  to  apply  a  linear  smoothing  formula,  such  as 
that  introduced  by  Longuet-Higgins  and  Cokelet  [Ref.  37]  in  their 
work  on  nonlinear  water  waves.  The  subsequent  applications,  includ- 
ing the  one  discussed  herein,  have  shown  that  chaotic  motion  sets  in 
sooner  or  later  regardless  of  the  smoothing.  The  repositioning  tech- 
nique introduced  by  Fink  and  Soh  [Ref.  40]  and  used  subsequently  by 
Sarpkaya  and  Shoaff  [Ref.  41]  removes  the  most  unstable  mode  and 
reduces  the  growth  of  the  higher  modes  of  Helmholtz  instability. 
However,  it  does  not  prevent  the  growth  of  spurious  waves  along  the 
vortex  sheet.  The  disadvantage  of  the  smoothing,  either  through  the 
use  of  a  numerical  filter  or  through  the  repositioning  of  the  vortices,  is 
that  it  is  not  clear  in  general  what  the  relationship  is  between  the 
results  achieved  and  the  unknown  exact  solution.  The  problem  may- 
not  possess  a  solution  for  all  time,  and  in  this  case  the  use  of  smooth- 
ing could  yield  an  acceptable-looking  solution  where  none  in  fact 
exists.  Alternatively,  the  solution  arrived  at  through  smoothing  may 
not  be  even  close  to  the  exact  solution. 
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B.     APPLICATION  OF  THE  DISTRIBUTED  VORTEX  METHOD 

The  liquid  surface  was  represented  by  a  number  of  point  vortices. 
Their  positions  were  sjonmetrical  with  respect  to  the  x  axis,  situated 
on  the  free  surface.  However,  the  strength  of  a  vortex  at  (x,  y,  nor- 
malized by  bo)  was  opposite  to  that  of  a  vortex  (image  vortex)  at  (-x,  y). 
From  a  mathematical  point  of  view,  one  would  like  to  have  the  vortices 
extend  from  -<»  to  +«»  and  the  two  trailing  vortices  originate  at  (1/2, 
-«o)  and  at  -1/2,  -^).  This  is  impossible  from  a  numerical  point  of 
view.  Observations  have  shown  that  the  free  surface  rises  in  a  rela- 
tively small  region  directly  above  the  trailing  vortices.  The  remainder 
of  the  free  surface  remains  undisturbed.  These  observations  and 
several  sample  calculations  led  to  the  conclusion  that  the  free  surface 
can  be  restricted  to  a  region  extending  from  x  =  -10  to  x  =  +10. 
Furthermore,  the  trailing  vortices  are  assumed  to  originate  at  a  depth 
of  y  =  -5. 

The  complex  function  representing  the  two  trailing  vortices 
below  the  free  surface  and  the  vortices  on  the  free  surface  is  given  by 

^Fq  iFo         .        _   .      iFm  ,      ,  > 

w  =  -  2^  In  (z  -  Zo)  +  2^  In  (z  +  Zq)  +  "2^  In  (z  -  zi) 

in  which  the  strengths  of  the  free-surface  vortices  are  normalized  by 
the  strength  of  the  trailing  vortex.  All  coordinates  are  normalized  by 
bo-  Then  the  velocities,  normalized  by  Vq  =  (ro/27ibo)  anywhere  in  the 
fluid  medium  is  given  by 
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_•         dw      _  ipo        1  ipo        1  irm        1  iFm        1 

^     ^^  -  dt  -      271  (z_z  J  ■"  2n  (z+ z  J  -^  271    (z-  zi)       27u    (z+ z^)     ^^^^ 


The  normalized  boundary  condition  at  the  free  surface  is  given  by 
D0m      Qm      Tim       ^ 


Thus,  at  least  theoretically  speaking,  one  can  calculate  the  poten- 
tial function  0  from  the  complex  velocity  potential  w  (the  real  part  of 
w),  the  velocity  of  the  vortices  from  equation  (35),  the  elevation  of  the 
free  surface  from  equation  (36)  for  a  given  Froude  number  Fv  (Fv  = 


Vo/Vgbo).  and  trace  the  evolution  of  the  free  surface  as  a  function  of 
time.  This  relatively  simple-sounding  procedure  is  anything  but  sim- 
ple, primarily  because  of  the  fact  that  the  numerical  instabilities  even- 
tually lead  to  large-scale  instabilities,  as  noted  in  the  introduction  to 
this  section.  The  use  of  various  smoothing  techniques  was  a  viable 
option  and,  in  fact,  was  tried  at  various  stages  of  the  calculations.  The 
results  have  shown  that  the  growth  of  the  instabilities  increases  with 
decreasing  Froude  number.  Neither  the  use  of  smoothing  techniques 
(e.g.,  the  Longuet-Higgins  technique)  nor  the  use  of  vortex  sheets, 
vice  discrete  vortices,  can  eliminate  the  eventual  development  of  a 
chaotic  behavior  on  the  free  surface.  It  is  because  of  this  reason  that 
the  use  of  smoothing  was  ruled  out  and  the  calculations  were 
restricted  to  a  relatively  high  Froude  number  (Fv  =  1.125). 

The  specific  details  of  the  numerical  calculations  are  as  follows: 
(1)  assign  the  position  of  the  vortices:  (2)  fmd  the  strength  of  the  vor- 
tices and  the  velocity  potential  assuming  the  free  surface  to  be  rigid: 
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(3)  advance  the  position  of  the  trailing  vortices  for  a  time  interval  t 
(0.01  in  the  calculations);  (4)  recalculate  the  velocity  of  the  surface 
vortices  and  advance  them  forward  in  time  for  a  time  interval  t;  (5) 
calculate  the  strength  of  the  free  surface  vortices  by  iteration  until  the 
free  surface  condition  Is  satisfied;  (6)  calculate  the  velocity  of  the  sur- 
face and  trailing  vortices  and  advance  their  positions  for  a  time  inter- 
val t;  and  (7)  repeat  the  calculations  by  returning  to  step  (5). 

The  procedure  described  above  is  relatively  simple  and  does  not 
require  excessive  computer  time  (about  20  minutes  on  IBM  3033). 
However,  the  free  surface  eventually  does  become  chaotic.  Some  of 
the  instabilities  can  be  alleviated  without  smoothing  through  a  judi- 
cious selection  of  the  initial  position  of  the  surface  vortices.  Numerous 
calculations  have  shown  that  the  vortices  near  the  y  axis  begin  to  slide 
sideways  as  the  free  surface  (or  the  vortex  sheet)  stretches.  Conse- 
quently, the  thinly  populated  regions  of  the  sheet  do  not  yield  a 
smooth  free  surface  and  the  flow  begins  to  leak  between  the  vortices. 
Also,  the  regions  where  the  free  surface  is  depressed  (where  the  scars 
develop)  become  overpopulated  with  vortices,  leading  to  the  growth  of 
short -wavelength  Helmholtz  instability.  This  problem  can  easily  be 
alleviated  by  packing  the  vortices  more  densely  near  the  y  axis  at  the 
start  of  the  calculations  so  that  the  entire  surface  becomes  more  or 
less  uniformly  represented  at  later  times.  This  simple  technique  has 
been  used  in  the  results  presented  herein.  An  exponential  function 
was  employed  to  assign  the  initial  vortex  spacings  so  that  near  the 
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y  axis  the  vortices  were  separated  by  a  distance  of  about  0.03  and  by  a 
distance  of  about  0.25  towards  the  end  of  the  vortex  sheet. 

C.     RESULTS  OF  THE  NUMERICAL  CALCULATIONS 

The  evolution  of  the  free  surface  with  time  is  shown  in  Figures  2. 1 
through  2.27  for  FV  =  1.125.  The  normalized  time  in  these  plots  is 
given  by  T  =  (Vq  t/bo  -  do /bo)  where  do  is  the  initial  depth  of  the 
trailing  vortices.  The  time  T  =  0  corresponds  to  that  at  which  the 
trailing  vortices  would  have  reached  the  undisturbed  free  surface  had 
they  continued  to  move  with  their  initial  mutual  induction  velocity. 
The  use  of  other  normalized  times  is  not  suitable  since  they  tend  to 
depend  on  the  initial  position  of  the  trailing  vortices. 

Figures  2.1  through  2.27  show  that  the  free  surface  directly  above 
the  vortices  rises  rapidly  while  the  adjacent  portions  of  the  surface 
depress  and  form  two  strong  scars.  Unlike  the  rigid  surface  case, 
where  the  trailing  vortices  continue  to  move  towards  right  and  left  at  a 
depth  of  about  y  =  -0.5,  the  trailing  vortices  below  the  deforming  sur- 
face almost  come  to  rest  at  a  point  slightly  above  the  free  surface. 
Furthermore,  their  initial  spacing  remains  nearly  constant.  This  sug- 
gests that  the  Kelvin  oval  formed  by  the  trailing  vortices  pushes  the 
free  surface  up  as  if  it  were  rigid  during  its  upward  migration.  Evi- 
dently, this  finding  is  based  on  the  assumption  that  the  trailing  vor- 
tices are  point  vortices  and  are  not  subjected  to  viscous  and  turbulent 
diffusion.  In  reality,  the  vortices  quickly  become  turbulent,  their  core 
radius  increases,  and  the  vorticity  is  diffused  over  an  ever-increasing 
area  with  the  passage  of  time.   The  amount  of  diffusion  depends  on  the 
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initial  position  of  the  trailing  vortices.  Consequently,  the  trailing  vor- 
tices emanating  from  large  depths  diffuse  over  a  larger  area  relative  to 
those  emanating  from  smaller  depths.  There  is  at  present  no  suitable 
mathematical  means  to  deal  with  the  turbulent  diffusion  of  such  vor- 
tices. The  best  one  can  do  is  to  generate  the  trailing  vortices  at 
depths  sufficiently  close  to  the  free  surface  to  prevent  excessive  dif- 
fusion and  yet  sufficiently  far  so  that  the  free  surface  remains 
undeformed  at  the  start  of  the  motion.  Extensive  calculations  and 
experimental  observations  have  shown  that  the  free  surface  does  not 
begin  to  deform  until  the  trailing  vortices  reach  a  depth  of  about  y  = 
-1.  Thus,  it  is  perfectly  safe  to  place  the  trailing  vortices  at  y  =  -3  at 
the  start  of  the  calculations. 
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V.  EXPERIMENTS 

A.  EXPERIMENTAL  APPARATUS  AND  PROCEDURES 

Experiments  were  conducted  in  a  water  basin.  It  consisted  of  a 
12-foot-long,  3-foot-wide  tank  with  aluminum  bottom  and  plexiglass 
walls  (see  Figure  3).  Additional  equipment  consisted  of  plumbing  for 
the  filling  and  emptying  of  the  tank,  a  coUimated  light  source,  and  the 
dye  system  for  flow  visualization. 

The  two-dimensional  trailing  vortices  were  generated  through 
the  use  of  two  counter-rotating  plates  (see  Figure  3).  The  mechanism 
rotating  the  plates  was  located  at  the  bottom  of  the  tank  and  below  a 
plexiglass  plate  spanning  the  entire  tank.  In  other  words,  the  motion 
of  the  driving  system  did  not  disturb  the  flow  above  the  plexiglass. 
The  mechanism  was  actuated  by  releasing  a  weight  attached  to  the 
common  driving  shaft.  Fluorescent  dye  was  introduced  slowly  into  the 
region  between  the  plates  through  the  use  of  two  holes  connected  to 
two  dye  reservoirs.  The  reason  for  the  use  of  two  reservoirs  was  that 
different  colors  of  dye  can  be  used  to  visualize  each  trailing  vortex. 

B.  PROCEDURES 

Experiments  were  initiated  by  filling  the  tank  to  a  suitable  level, 
removing  any  air  bubbles,  bringing  the  plates  to  their  full  open  posi- 
tions, waiting  for  a  sufficient  period  of  time  for  the  fluid  to  come  to 
rest,  introducing  the  neutrally  buoyant  dye,  actuating  the  light  source 
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and  the  video  system,  and  initiating  the  rotation  of  the  plates  by  sud- 
denly releasing  the  load  attached  to  the  driving  shaft.  The  plates 
rotated  smoothly  and  then  came  to  rest  on  the  plexiglass.  In  other 
words,  the  plates  generating  the  vortices  literally  "disappeared." 
Thus,  no  other  vortices  were  generated.  It  is  a  well-known  fact  that 
this  is  not  the  case  v^th  piston-generated  vortices.  When  the  piston 
stops,  two  additional  vortices  (from  a  two-dimensional  piston)  or 
another  vortex  ring  (from  an  axisjnnmetric  piston)  are  generated.  The 
mechanism  used  in  the  present  investigation  effectively  prevents  the 
generation  of  secondary  vortices  by  simply  disappearing. 

The  trailing  vortices  rise  under  their  mutual  induction  velocity, 
quickly  give  rise  to  a  Kelvin  oval,  and  smoothly  approach  the  free  sur- 
face. When  the  vortices  are  at  a  distance  of  about  1.0  from  the  free 
surface,  the  free  surface  begins  to  rise  and  forms  a  smooth  hump,  bor- 
dered by  two  scars.  For  the  Froude  numbers  achievable  (maximum 
0.35),  the  trailing  vortices  turn  quickly  to  right  and  left  (as  if  they 
were  approaching  a  rigid  surface),  the  scar  front  moves  in  the  respec- 
tive directions  ahead  of  the  vortex,  and  the  hump  between  the  vortices 
begins  to  recede.  During  the  later  stages  of  the  motion,  the  scars  are 
slaved  to  the  vortices  and  continue  to  move  ahead  and  in  the  direction 
of  the  vortices. 

Figures  4.1  through  4.12  show  the  evolution  of  the  free  surface  at 
various  times  T  for  a  Froude  number  of  Fv  =  0.35.  The  vortex  centers 
are  clearly  visible  in  these  pictures  because  of  the  additional  assistance 
provided  by  the  dissolved  air  bubbles  to  the  flow  visualization  efforts. 
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VI.    DISCUSSION  OF  RESULTS 

The  computer  code  based  on  the  boundary  element  methods 
through  the  use  of  the  discrete  vortices  has  been  used  to  carry  out 
calculations  for  various  Froude  numbers  in  order  to  predict  the  evolu- 
tion of  the  free  surface  deformation.  The  t^npe  of  the  vortex  distribu- 
tion, time  increment,  and  the  iteration  techniques  have  been  varied 
within  limits  to  explore  the  differences  in  the  motion  of  the  free  sur- 
face. The  results  have  shown  that  the  calculations  are  fairly  stable  at 
relatively  high  Froude  numbers.  However,  at  relatively  small  Froude 
numbers,  the  Kelvin-Helmholtz  instability  sets  in  quickly  and  the  free 
surface  exhibits  highly  irregular  shapes.  It  appears  that  either  the 
integration  procedures  have  to  be  refined  or  more  sophisticated  vor- 
tex sheets  have  to  be  used  to  calculate  the  small  deformation  of  the 
free  surface  at  small  Froude  numbers. 

The  results  obtained  with  a  Froude  number  of  Fv  =  1.125  are 
showTi  in  Figures  2.1  through  2.27.  The  most  striking  feature  of  these 
results  is  that  the  free  surface  rises  rapidly  to  a  height  of  about  1.25 
above  the  mean  water  level  and  captures  the  vortex  pair.  It  has  not 
been  possible  to  carr^^  out  the  calculations  to  larger  times.  The  fact 
should  be  kept  in  mind  that  the  point  vortices  used  in  the  numerical 
calculations  may  have  ver^^  little  resemblance  to  the  real  vortices  by 
the  time  they  rise  to  the  mean  water  level.     In  fact,  the  experiments 
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show  that  the  vortices  become  quickly  turbulent.  There  is  at  present 
no  possibility  of  incorporating  into  a  numerical  analysis  the  motion  of  a 
turbulent  vortex. 

Experiments  were  carried  out  at  a  maximum  Froude  number  of 
0.35.  The  evolution  of  the  free  surface  is  showTi  in  Figures  4.1  through 
4.12.  The  shape  of  the  free  surface  at  the  time  of  maximum  rise  is 
plotted  in  Figure  5.  Figures  4  and  5  show  that  the  free  surface  rises  to 
a  maximum  height  of  about  0.25.  develops  two  strong  scars,  and  then 
subsides  as  the  vortices  begin  to  move  parallel  to  the  free  surface.  The 
scar  front  is  slightly  ahead  of  the  vortex  center  and  reduces  to  a  small 
depression  at  larger  times.  The  results  presented  in  Figures  4  and  5 
should  form  the  basis  of  comparison  for  future  numerical  efforts  at  the 
corresponding  Froude  numbers.  Additional  calculations  are  underway 
with  more  sophisticated  vortex  sheets.  The  results  presented  herein 
are  extremely  encouraging  and  are  expected  to  lead  to  a  better  under- 
standing of  this  extremely  complex  and  challenging  phenomenon. 
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VII.    CONCLUSIONS 

The    investigation    reported    herein    warranted    the    following 
conclusions: 

1.  The  basic  equations  governing  the  motion  of  vortices  in  homo- 
geneous and  density- stratified  media  have  been  developed  and 
expressed  in  terms  of  dj^amic  and  buoyant  scaling  laws; 

2.  A  novel  experimental  technique  has  been  devised  to  generate  a 
nearly  two-dimensional  vortex  pair  rising  toward  a  free  surface; 

3.  Migration  of  a  vortex  pair  toward  a  free  surface  gives  rise  to  a 
bulge  and  two  scars; 

4.  The  maximum  rise  at  the  free  surface  occurs  when  the^  vortex 
pair  is  at  a  depth  of  about  0.5; 

5.  In  the  range  of  Froude  numbers  encountered  in  the  experi- 
ments, the  rise  of  the  free  surface  is  always  followed  by  a  fall  as 
the  vortices  begin  to  move  parallel  to  the  free  surface  at  large 
times; 

6.  The  numerical  analysis  of  the  free  surface  deformation  through 
the  use  of  discrete  vortices  leads  to  instabilities  at  low  Froude 
numbers.  These  instabilities  could  have  been  removed  with  a 
suitable  numerical  filter  or  smoothing  technique.  However,  the 
use  of  a  relatively  subjective  smoothing  procedure  has  been  ruled 
out.  The  calculations  at  a  Froude  number  of  1.125  have  yielded 
fairly  stable  solutions  and  provided  the  basis  for  comparison  with 
future  experiments. 
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APPENDIX:    FIGURES 
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Figure  1.         Computational  Domain.  Governing  Equations, 
and  Boundary  Conditions 
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Figure  2.1      Velocity  Field  for  T*  =  -2.40 
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Figure  2.2      Velocity  Field  for  T*  =  -2.30 
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Figure  2.3      Velocity  Field  for  T*  =  -2.20 
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Figure  2.4 


Velocity  Field  for  T*  =  -2.10 
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Figure  2.5      Velocity  Field  for  T*  =  -2.00 
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Figure  2.6      Velocity  Field  for  T*  =  -1.90 
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Figure  2.7      Velocity  Field  for  T*  =  -1.80 
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Figure  2.8      Velocity  Field  for  T*  =  -1.70 
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Figure   2.9      Velocity  Field  for  T*  =  -1.60 
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Figure    2.10   Velocity  Field  for  T*  =  -1.50 
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Figure   2.11    Velocity  Field  for  T*  =  -1.40 
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Figure  2.12   Velocity  Field  for  T*  =  -1.30 
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Figure   2.13   Velocity  Field  for  T*  =  -1.20 
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Figure  2.14   Velocity  Field  for  T*  =  -1.10 
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Figure  2.15   Velocity  Field  for  T'  =  -1.00 
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Figure  2.16   Velocity  Field  for  T*  =  -0.90 
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Figure  2.17   Velocity  Field  for  T*  =  -0.80 
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Figure  2.18   Velocity  Field  for  T*  =  -0.70 
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Figure  2.19   Velocity  Field  for  T*  =  -0.60 
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Figure  2.20   Velocity  Field  for  T*  =  -0.50 
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Figure  2.21    Velocity  Field  for  T*  =  -0.40 
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Figure  2.22   Velocity  Field  for  T*  =  -0.30 
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Figure  2.23   Velocity  Field  for  T*  =  -0.20 
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Figure  2.24   Velocity  Field  for  T*  =  -0.10 
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Figure  2.25  Velocity  Field  for  T*  =  0.00 
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Figure  2.26    VelociU^  Field  forT*  =  0.10 
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Fmure  2.27   Velocin^  Field  for T*  =  0.20 
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